#include <stdio.h>
#include <stdlib.h>
#include <sys/time.h>

#include <omp.h>

double getPi(int iterations, int threads);

int main(int argc,char **argv)
{
    int iterations = 0;
    double time[3];
    double result = 0.0;
	int nThreads = 0;
	const int measure_iterations = 1000;

	/* set options */
    if (argc == 3)
	{
		iterations = atoi(argv[1]) >= 1 ? atoi(argv[1]) : 100;
		nThreads = atoi(argv[2]) >= 1 ? atoi(argv[2]) : 2;
    }
	else
	{
		printf("Usage: %s <iterations> <threads>\n", argv[0]);
        return -1;
	}

	/* calculate pi */
	time[0] = omp_get_wtime();
	for (int i=0; i<measure_iterations; i++)
	{
		result = getPi(iterations, nThreads);
	}
	time[1] = omp_get_wtime();

	/* print result */
	time[2] = (time[1]-time[0]) / measure_iterations;
    printf("%i %i %lf %19.18lf\n", nThreads, iterations, time[2], result);

    return 0;
}

double getPi(int iterations, int threads) // mhmmmm... pie...
{
	int i;
	double res_thread, res_global;
	double deltaX;
	double xi;

	deltaX = 1.0 / (double)iterations;

	res_global = 0.0;
	omp_set_num_threads(threads);
	#pragma omp parallel private(i, res_thread) shared(res_global)
	{
		res_thread = 0.0;
		#pragma omp for
		for (i=0; i<iterations; i++)
		{
			xi = (double)i / (double)iterations;
			res_thread += (4.0/(1.0 + xi*xi))*deltaX;
		}
		//printf("%2d: rt=%lf\n", omp_get_thread_num(), res_thread);
		#pragma omp atomic
		res_global += res_thread;
	}

	return res_global;
}

